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Abstract: Within a finite element context, the embedded strong dis- 
continuity approach for strain localization has been used extensively to 
model localized deformation and fracture in geonraterials. As a fracture 
propagates, any changes in orientation will inhibit sliding and force the 
surfaces, in some locations, to open. Some previous models feature only 
a single sliding degree of freedom. As the surfaces slip in such models, an 
artificial hardening occurs, creating a geometric locking effect. To this 
end, we implement a formulation, which, in addition to sliding, possesses 
an opening degree of freedom. We develop a constitutive model that al- 
lows for coupled opening and shearing displacement in tension, as well 
as pure sliding in compression. This paper compares the single degree of 
freedom formulation with the model containing both degrees of freedom. 

We show that the locking effect is alleviated. In addition, a method for 
increasing the robustness of softening problems is featured. 

Keywords and phrases: strong discontinuity, enhanced strain finite 
element, bifurcation, localized deformation, simulation robustness. 



1. Introduction 



The numerical modeling of geonraterials, especially in the context of finite 
element analysis, has a wide and varied base of research behind it. As part 
of these developments over the past two decades, embedded discontinuities 
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have been investigated for use in numerical models which seek to capture 
failure due to localized softening. The concept of a so-called weak disconti- 
nuity was developed in order to capture narrow shear bands of deformation 
within materials [3, 38]. Shear bands were first described in metals by Tresca 
[39]. Later, this nomenclature was adopted by the geomechanics community 
in order to describe similar phenomenon seen in soil and rock-like materials. 
Subsequently, numerical models were developed in order to replicate this be- 
havior in simulations. Precipitated from this line of research was the strong 
discontinuity approximation. When the width of the discontinuity is small 
compared to that of the overall structure the weak discontinuity may be sim- 
plified by reducing it to zero width. Simo and co-workers showed in seminal 
reports [36, 37] that the post-localization softening modulus of the material 
as well as the plastic multiplier necessarily take on a distributional represen- 
tation. In addition to this representation, they introduced a re-parametrized 
displacement and strain field which is suited for capturing both measures of 
deformation within a finite element context. 

The model presented in this paper is an extension of the Mohr-Coulomb slip 
model for geological materials in [5, 10, 13, 31, 32]. The main contribution 
of this paper is the introduction of a post-localization constitutive model 
allowing for slip in directions both tangential and normal to a given disconti- 
nuity band. The impetus for this is two-fold. In reality, as geomaterials begin 
to localize and subsequently fail in either the compressive or tensile regime, 
material may not only slide but also, noticeable movement in the direction 
normal to the failure surface is observed. Even for far-held compression, the 
development of a model which takes into account the interaction of both 
compressive and tensile traction necessary. 

Carol et al. [8] developed a hyperbolic yield surface in order to capture the 
effect of tension and a linear one for the purely compressive regime. This 
model was later adopted by Galvez, et al. [15] for the purposes of simulating 
material failure in masonry structures. The model in this paper is similar 
to that of Carol and Galvez in that it captures the aforementioned defor- 
mation regimes. Camacho and Ortiz [7] developed a model which features 
an elliptical interaction between the shear and normal tractions on a given 
deformation band. In this model, a material parameter is employed to con- 
trol the contribution of the shear traction on the surface, a feature which 
is later used in [40]. In a similar a vein, we develop a model which features 
elliptical interactions between both the tractions and displacements (normal 
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and tangential to the surface cf. [41]), based on the material’s strain energy 
release rate. 

Secondly, in [13] simulations yielded results which revealed that a model 
only incorporating a single, sliding degree of freedom gave rise to geometric 
locking which was mesh-dependent. Specifically, this effect arose due to the 
model’s inability to accurately capture opening displacements. From this ob- 
servation, the conjecture was that remedies such as extra degrees of freedom 
[23, 26] would alleviate geometric locking. The authors in this paper took 
the latter approach by developing a constitutive model which features both 
opening and sliding degrees of freedom on the discontinuity surface. When 
the band is in compression, it is only subject to deformation tangential to its 
surface. Hence, friction and cohesion forces impede sliding on the localization 
band. However, when tension is present the band experiences a combination 
of opening normal to the surface in addition to the aforementioned tangen- 
tial sliding, albeit without friction. We show that such a model presents a 
more realistic depiction of the mechanics of softening behavior within soil- 
based materials. Through several benchmark examples, this paper shows that 
the currently proposed constitutive model alleviates the spurious geometric 
locking effect seen in [13]. 

In this particular study, we focus on a damage-like formulation which as- 
sumes elastic unloading toward the origin and the tradition of cohesive zone 
formulations. For simplicity, we will consider a constant friction coefficient. 
However, this is not a requirement. By ramping up the friction coefficient as 
the cohesion degrades, we can recreate the simplified slip weakening formu- 
lation as described in Borja and Foster [5]. Variable friction as detailed in 
[12] and [13] can also be implemented. 

Softening problems within the strong discontinuity framework are well known 
to have robustness and stability issues [4, 20, 27]. To this end, we incorporate 
a combined implicit-explicit integration method (Impl-Ex), modified from a 
scheme devised by Oliver and co-workers [28], which renders the local stiffness 
matrices associated with softened elements positive definite. We show that 
this allows problems with complex geometries to converge, whereas using a 
fully implicit integration scheme does not. 

The remainder of the paper is structured as follows: Section 2 gives a general 
overview of strong discontinuity kinematics in the context of small strain. 
Additionally, in this section, the localization criterion, which indicates the 
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onset of the formation of a strong discontinuity, is given a brief review. Section 
3 briefly discusses some fundamental aspects of an enhanced strain finite 
element. The constitutive model for slip weakening within a localized element 
is described in Section 4. The numerical implementation of the model is given 
in Section 5. Section 6 shows the derivation of the stiffness formulation. 
Numerical examples and closing remarks conclude the paper. 

A summary of the symbolic notation used throughout this paper is noted 
here. Boldface quantities denote vectors and tensors, e.g., (a) y . = a lt] . The 
symmetric dyadic product in index notation is (a <8) b)®. = | {a>%bj + ajbi ), 
where the superscript “s” denotes a symmetric second-order tensor. The dot 
product is expressed as a - b — a t b t , where Einstein’s summation convention 
is used. The inner product symbol is used for contraction of two second- 
order tensors as a : b = ayby or for the contraction of a fourth-order tensor 
and a second order tensor as c : a = c y/ y a^i . The symmetric gradient operator 
acting on a vector is written as (V s a)- = \ ( a,ij + Oj t i). 



2. Kinematics 

2.1. Continuum Equations 

Here we describe the kinematical relationship between continuous and dis- 
continuous displacements within an arbitrary continuous volume (Figure 1). 
Note that for the purposes of this study, our analysis is relegated to the 
realm of infinitesimal strains. The continuous displacements are distributed 
throughout the volume while the discontinuous displacements are contained 
on the discontinuity surface S. The total displacement held in the body is 
given by the form: 



u := u+lu}H s (x) (2.1) 

where u is the vector of continuous displacements and [it] designates the 
vector of displacements along the strong discontinuity. This vector takes the 
form: 



[«] = < = + c»« 



( 2 . 2 ) 
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Fig 1: Arbitrary volume with embedded strong discontinuity surface 

Where l and n are unit vectors which are tangential and normal to the dis- 
continuity surface, respectively. 

Hs ( x ) is the Heaviside function across the discontinuity surface S defined 
by the conditions 



The strain tensor is found by taking the symmetric part of the gradient of 
the displacement vector: 



where VH$ = nSs ■ The physical meaning of the last term is that the strain 
is unbounded on the discontinuity surface S. 

For the purposes of this study, the jump magnitude is considered to be spa- 
tially invariant (piece-wise constant across the elements). Thus the second 
term is omitted from the equation, leaving the final form of 




x e Q + 
x e 



( 2 . 3 ) 



e := V'u = V*u + V" [u] H s + (M ® n)‘ S s 



( 2 . 4 ) 



e V°u — V s u + ([«] ® n) s S s 



( 2 . 5 ) 
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2. 2. Localization Condition 

In this section, the procedure that is used to detect the onset of strain local- 
ization is outlined. Hill [17] investigated the physics of wave propagation in 
solids in order to determine the onset of inelastic behavior therein. Rudnicki 
and Rice, in [35], built off of this work introducing a mathematical frame- 
work for detecting shear band localization. This formulation has since been 
extended to the case of displacement discontinuities as well as strain disconti- 
nuities ([6, 30] among others). Since then, bifurcation theory has been widely 
used to determine the onset of softening behavior of various geomaterials. 

To begin, we posit that at any given time, the traction rate along the discon- 
tinuity surface must be continuous. Using a general non-associative plasticity 
model, i.e. , where F is an arbitrary yield function and G is a plastic potential, 
the traction rate is written as 



t = & ■ n = n ■ c e 




(2.6) 



Where c e is the elastic modulus tensor and A is the rate of the plastic mul- 
tiplier. 



For a discontinuity in the displacement field, the plastic multiplier must have 
a distributional form. From [6] A = A $83 

where \ 5 = : c e : ([ft] ® n) S and X = f£ : c e : §§ 

Substituting this into (2.7) yields 

t = n ■ c e : V s R y + c ep : (fwj <g> n) S 5 S (2.7) 

^ V ^ ✓ 

N/* 

1 i s 

where the term c ep = c e — A c e : ^ <g) J2A : c e is the elastic-perfectly plastic 
tangent modulus. 
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In order for the traction to be bounded, the discontinuous part must vanish, 
hence 

t 5 = n ■ c ep : ([it] <g> n) s = 0 

Noting that [it] = ( m , this expression can be re-written as 

(n • c ep ■ n) ■ m = A ■ m = 0 (2.8) 

where A is known as the elastic-perfectly plastic acoustic tensor. 

This quantity is zero only when A becomes singular, that is, when det ( A ) = 
0 . 

Eq. 2.8 can be solved iteratively, as detailed in [29] and [42], in order to 
determine the direction of the jump, m, along the discontinuity. 



3. Post-localization Constitutive Model 
3.1. Yield, Criteria 

Once localization has been detected and the orientation of the critical surface 
determined, the mechanics of the softening behavior on the surface must be 
defined. There are several ways to approach this model. Oliver [25] has shown 
that the continuum model can be used to induce a localized model. The great 
advantages of this approach are that a smooth and consistent transition 
from continuum to localized response is assured, and the material response 
is consistent with experimental data. However, it may also be argued that the 
mechanics of the localized material differ from the bulk continuum response 
since this is a separate mechanism. In this case, a separate constitutive model 
must be introduced for the localized behavior. The form of that model is the 
subject of this section. 

The initial cohesion on the surface is determined to be consistent with the 
stress state at localization. The yield strength in pure tension may be differ- 
ent, and thus is assumed to differ by a constant factor of a a . The constant, 
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a a can be thought of as the reciprocal of the parameter /3, which is described 
as a shear stress factor in Camacho and Ortiz [7]. Similarly, Alfaiate et al. 
[1] use a parameter which controls the contribution of the shear jump on the 
deformation band. 

The yield function that we incorporate takes two forms. When compressive 
normal traction is present on the band, the yield function takes the form of 
a Mohr-Coulomb law (Eq. 3.1). In the case of tensile normal traction, the 
interaction between the shear and normal traction is assumed to be elliptical 
(Eq. 3.2). 



^compression 0 C • sign (<^ s ) | f (3-1) 

^tension = 0 = \J (jf + (a ff {(j)f ~ C (3.2) 

where a = n ■ cr ■ n is the normal traction, r = l ■ cr ■ n the shear traction 
on the deformation band. ( s is the amount of slip in the direction parallel to 
the deformation band and c is the cohesion. (•) are the Macaulay brackets 
[{x) = (x + |x|) / 2], signifying the positive part of the quantity. Hence in 
compression, the normal traction only influences the yield function in the 
frictional traction / = n(—a). For this study, we use a static coefficient of 
friction, though a variable coefficient, as in [5, 13], can be used. 



3.2. Cohesion Softening Formulation 

The cohesive formulation we will use takes the form of a traction- displacement 
relationship, common in the cohesive zone formulations such as [7, 11, 18]. 
This formulation follows the spirit of localized damage mechanics, so that 
the interface softens with increased displacement, but unloads elastically. 
Reloading is also elastic until the previous displacement level is reached, 
where damage continues (Figure 3). For our particular formulation, we use a 
linear form given by the equation 




C eq 

c 7 



C = Co 



(3.3) 
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Fig 2: Normal and shear traction interaction on yield surface 



where c is exactly the cohesion along the surface in standard Mohr-Coulomb 
criteria, cq is the initial cohesion at the time of localization, ( eq is a weighted 
magnitude of the slip on the band and (* is the slip distance until the cohesion 
completely degrades to zero. 

In the literature, several other cohesion softening laws, particularly suited for 
brittle materials such as concrete, are bilinear [16, 44], exponential [9, 43], and 
power- law based [33]. For geomaterials however, a linear model for cohesion 
degradation seems to fit experimental data well [19, 34], 

The initial cohesion is determined to be consistent with the bulk stress state 
at the moment of localization. In addition to linear softening we assume 
elastic unloading and reloading with respect to the equivalent traction and 
displacement. 

We define a parameter ( c = (* (co — c) /cq, which is the maximum equivalent 
slip magnitude observed on the slip surface up to the current time. The 
parameter k s is defined as the slope of the unloading-reloading curve at a 
given value of ( c . Therefore, k s = c/( c . On this portion of the curve, we also 
assume that the cohesion is frozen at a value of ( c , thus the explicit form of 
k s is given by 



D. A. Weed et al./A Novel Post-localization Model 



10 




Fig 3: Cohesion softening law displaying different unloading-reloading curves 
at arbitrary time steps n and n+1 



*■ = 

= *Gr?) 

Further softening on the band can only occur when ( eq > ( c . 

Remark: The assumption that the cohesion is constant on the unloading- 
reloading portion of the graph is a simplification but also one which is 
grounded in reality. One could imagine that, under very small displacement 
values, as the crack moves back and forth, some material will remain intact 
although it will stretch or elongate. It is quite possible for a cohesive mate- 
rial to undergo this process. For very small displacements, the material will 
deform but the overall cohesive force of the soil will remain intact. 



(3.4) 
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3.3. Tensile Regime 

In tension, the friction force is absent, therefore according to Eq. 3.2 the 
equivalent stress on the band is 



Oeq — c— \/t 2 + [0c o a) 



(3.5) 



We relate the stiffness parameter, k s , to the traction on the band 



r 2 + (a a a) 2 = c 



cr, 



eq 



I 1 C e <7 

“ 1 C 



C eq A Ceq 



C0 V (* ) Ce, 

“"“(we) 0, 

= eq 



(3.6) 



where ( eq = \J ( s 2 + («cCn) 2 - The material parameter, a^, signihes the differ- 
ing weights of the opening and shear displacements. 

We now postulate two balance laws which relate the tractions and displace- 
ments both normal and tangential to the deformation band: 



(3.7) 

(3.8) 



Eq. 3.7 can also be written as 

o' = —k s (n 

QL(j 



OL(jG ks^C ) Cn 

t = k,C„ 



(3.9) 
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where k n represents the stiffness in the direction of the normal traction on 
the band. If we rewrite this stiffness parameter, solving for a a , we have 
ot a = Trk s . This is substituted into Eq. 3.5 giving 




Further, substituting the forms in the traction- displacement balance laws in 
Eqs. 3.7 and 3.8 yields 



& eq y (k s ( s ) T (oy^'sCn) 

= hy/Ca + («cC nf 

ksCeq 

which recovers Eq. 3.6, verifying the coupled nature of the tractions and 
displacements on the deformation band. 





Fig 4: For a deformation band experiencing tension, the slip degrees of free- 
dom are coupled with the tractions. 



3.3.1. Determining the parameters a a and a ^ 

These parameters are related to one another vis a vis the fracture energy in 
each of the respective fracture modes. The specific fracture energy for Mode 
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I (opening) and Mode II (sliding) (Figure 5) is simply the area under the 
respective curves (cf. [45]). Therefore 




Fig 5: The fracture energy for opening and sliding fractures 



The fracture energy is related to the stress intensity factor, K, through the 
equation 



G = ^(l-i/ 2 ) (3.12) 

which is valid for plane strain. 

Hence, the ratio of the corresponding stress intensity factors is 



K 



ii 



Kj 



G 

~G, 



ii 



2 



oiq ■ a a 



(3.13) 
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There are many solutions to Eq. 3.13, however in this study, we choose the 



K 



solutions in which = a a . In [2], the empirical value for the ratio — — varies 

Ki 

between 1.13, which corresponds to softer marble, and 2.19, for limestone. 
In this study, we simply choose = a a ~ 2.0. In a future study that will 
feature the combination of this model and the plasticity model detailed in 
[24], we choose a value closer to that of limestone (which is ~ 2.14), for that 
will be the particular material of interest. It is worth noting that authors in 
[7] derived = 1 / a a on energy considerations. However, this creates equal 
fracture energy in Mode I and II, which is not observed in many materials. 



3-4- Compressive Formulation 

The cohesive formulation is exactly the same as in the tensile case. However, 
the normal effects on a eq are neglected. In other words, under compression 
cr eq = t. In addition to this, the normal displacement on the band is assumed 
to be zero. Hence, r = k s ( s . 

However, in addition to cohesive tractions, there is a friction force, defined 
earlier as / = Under compression, the frictional force acts indepen- 

dently of the displacement, and hence is always active. The final form of the 
traction balance in this case may be written as 



\r-k s Cs\<f (3.14) 

where strictly less than implies no motion on the band, and equality allows 
for slip. 

For many quasi-brittle materials, there is actually noticeable dilation along 
shear bands and fracture interfaces due to the mismatch of rough surfaces, 
wear, and the resulting gouge material [22] . However, this often small and is 
beyond the current scope of this work. 
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4. Finite Element Implementation 
4-1. Re-parameterization of displacement field 

Here, we use the parametrized form of the displacement introduced by Simo 
and co-workers [37]. Within a standard Assumed Enhanced Strain (AES) 
framework, the displacement held is recast such that the continuous part 
conforms to the standard finite element formulation and the discontinuous 
part is considered as an additional enhancement to the displacement held. 
Hence, 



,h,conf ^jh,enh 


(4,1a) 


+ [[V]] Mg (X) 


(4.1b) 



where Mf ( x ) = H s — / h . Hence, 

u h = (u h + [[u h ]} f h ) + \[u h ]} (Hs - f h ) (4.1c) 



The term, f h , is an arbitrary smooth function which, on a given element with 
a discontinuity, meets the criteria 



f 1 if e n E 
( 0 if e n E fl~ 



(4.2) 



and e n is the element node number. For this study, f h takes the particular 
form of 



^en 

f h = Y J N A H s {x A ) 

A =i 



(4.3) 
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where n en is the set of active nodes within a given element, i.e. nodes that 
are in Q + . For each element, the conforming displacement held is composed 
of the continuous nodal displacements as well as the contribution of the 
displacements along the localization band. This part of the displacement 
“conforms” to the standard finite element shape functions, hence 



hen 

u h,conf = N A d A (4.4) 

A= 1 

The strain tensor is obtained by taking the symmetric part of the gradient 
of the displacement in eq. (4.1c) yielding, 



e h = Bd + [(C ® n) s 5 S - V/ ft ) *] (4.5) 



^ h,conf _|_ ^ h,enh 

In terms of recovering the finite element stress, it is convenient to denote the 
strain in terms of a regular or continuous part and a singular or jump part: 

e h = Bd- (C^V/Y+ ^nffe (4.6) 

regular singular 



as shown in [6], the finite element stress for localized elements is derived from 
the regular part of the strain, hence 

& h = c e : e h ’ re9 (4.7) 



5. Numerical Implementation of Slip Model 

The solution technique employed to solve the slip at the element level is a 
piecewise Newton-Raphson iteration. This particular approach is well-suited 
for the proposed model due to the fact that the combined opening-sliding 
model is accompanied by a nonlinear solution. An iterative scheme is also 
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a good choice because it has the potential to capture a broader range of 
constitutive models that may be the focus of future investigation. 

Unfortunately, it is not always possible a priori to determine whether the 
band is in compression or tension. This adds a non-smoothness to the equa- 
tions that, in some cases, creates solution difficulties. 

To this end, the slip parameter is solved using two distinct subroutines. One 
subroutine handles the case of positive normal traction on the band, thus 
solving for both opening and shear slip parameters. If the band is in com- 
pression (hence, zero normal slip on the band), a subroutine which solves 
solely for shear slip is utilized. If a change in the sign of the normal traction 
is detected within a given subroutine, then a new slip value is interpolated 
and used as the starting value in the appropriate subroutine. Doing this en- 
sures that the Newton-Raphson is settling on the correct solution, given the 
traction state on the band. 



The two subroutines differ in terms of which traction balance equations are 
utilized. Namely, when the band is in compression, only one equation is 
needed, which is the balance between the shear traction, friction and cohesive 
forces. In this case, the shear slip is solved for via the equation 

$ = 0 = \t - /csCsl - /• (5.1) 

Even if the N-R successfully yields a converged slip value, it must be ensured 
that the normal traction sign is negative (hence, the band is in compression) 
since it is possible for the band to go into a state of tension. This simply 
indicates that the slip on the band should instead be calculated using the 
combined opening-sliding formulation. In order to find the value of ( s at the 
onset of tension (hence, when cr = ( n = 0.0) a simple linear interpolation is 
performed using Eq.5.2. We then use these new values of ( n = 0.0 and ( s to 
initiate a N-R iteration in the opening-sliding subroutine (Eqs. 5.3 and 5.4). 
The three points used for interpolation are as follows: (1) (a®, £*), signify the 
initial normal traction and slip values at the beginning of the N-R iteration, 
(2) (cr = 0.0, Cs), the unknown value of ( s when the normal traction is zero 
and finally (3) (cry, £/) which are the spurious converged values of the positive 
normal traction and the shear slip. 
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C. 



c 



a 



c/ - c 






cr 



(5.2) 



When the band is undergoing tension, the friction force is absent, and only 
the normal and shear traction forces need to balanced. Thus the following 
set of balance equations is simultaneously solved for the normal and shear 
slip values: 



= 0 = a - k n C n (5.3) 

$ 2 = 0 = r - k s Cs (5.4) 

Again, after convergence, the normal traction sign must be checked. If it is 
negative (hence the band should be in compression) then the normal slip, ( n , 
will also be negative, which is of course a spurious value. In this case, a new 
value for ( s is calculated using the equation 



(5 ' 5) 

which is the value of the shear slip at the onset of compression (i.e., when 
a = ( n = 0.0). After this value is found, it is used as the initial start point 
for the N-R iteration in the sliding-only formulation (Eq. 5.1). 

The solution procedure for either formulation follows a standard Newton- 
Raphson approach given by 



A°k-\- 1 A>k 

>n+l >n+l 



38 <g + i) t 

aCi ) 



* (CO 



-i 



(5.6) 
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5 . 1 . Slip Algorithm 

The solution algorithm differs in form depending on whether the element is 
newly localized (in which case Box 1 is used) or if slip has already begun on 
the band (Box 2). 



Box 1: Slip algorithm for a newly localized element. 



Step 1 : Compute rr( l r +1 = er n + c e : Ae cor U 
Step 2: Check for yielding on the band: $ > 0? 

If no, band is inactive. Set 

<T n+l = cr 7i + l 

Cn+1 Cn 
and exit. 

Step 3: If d> > 0, proceed to solve for slip values on the band: 

If a > 0 (tension) solve for £ n and £ s using equations 

O' - k n ( n = 0 
T - k s ( s = 0 

After convergence, if cr < 0 then interpolate 

a new sliding value (C,) using Eq.5.5 and use 

this as a starting value for the sliding-only formulation. 

Else (cr < 0, hence compression) 

Solve for £ s using: |r - fc s £ s | - / = 0 

After convergence, if cr > 0 then interpolate 
a new sliding value using Eq.5.2 and use this 
as a starting value for the combined formulation. 

Step 4'- Update £ n+1 , k s n+ 1, and <x n+1 then exit. 
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Box 2: Slip algorithm for an clement which has pre-existing slip on the band. 



Step 1 : Compute tr% +1 =cr n +c e : Ae con f 

Step 2: Check the traction on the band and enter the appropriate 
subroutine. First, assume elastic unloading/reloading on the band, hence 
hold k s and k n as constant and proceed to solve the balance equation(s). 

If cr > 0 (tension) solve for £ n and £ s using equations 

cr ~ k n Cn = 0 

r - k s ( s = 0 

Else (cr < 0, hence compression) 

First, check for slip on the band: 

If | t - k s ( s \ < f 

No slip on the band due to friction, 
exit with trial stress. 

Else 

Solve for ( s using: |r — k s ( s \ — f = 0 

After slip value(s) calculated in elastic unloading/reloading 
phase, check: 

If Ceq < Cc 

It is confirmed that the band is either 
unloading or reloading. Therefore, update 
(Tn+l and set Cc,n+1 = Cc,n an( I ex if- 

Else 

The band is in the softening phase. 

Recalculate the balance equation(s) allowing k s and k n 
to vary as in Box 1. Additionally, interpolate as needed 
if the traction sign changes. 

Step 3: Update Cn+ii ^s,n+ i, fcn.n+i and <r n +i then exit. 
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6. Stiffness Matrix Formulation 

For the strong discontinuity approach, the resulting stiffness matrix from this 
formulation, assuming elastic unloading in the bulk material can be derived, 
following [13]. We begin with two sets of equations that must be solved: 
the standard balance of linear momentum (here taken to be quasi-static and 
small strain) and the traction balance on the localized deformation band. 



r e = B t . crdfi 



n% dn 






NHdT = 0 

$n+ 1 = 0 



( 6 . 1 ) 

( 6 . 2 ) 



Taking variations on these, we arrive at 



5r e = K e dd 5d e + K^C (6.3) 

5$ = K e Cd 5d e + K e (C 5C e (6.4) 

where K dd is the standard element stiffness matrix, and the others can be 
shown to have the following forms: 



K e dr = 


-i 


B‘,c V 


«c 




dC 






d® „ 


K e , d = 




= . c . B 


Qa 


~dd ~ 


da 


K tc = 


d $ 
9C 





(6.5) 

( 6 . 6 ) 

(6.7) 



The tensor B is the third-order, symmetric gradient of the nodal displace- 
ment interpolation functions, commonly referred to as the strain-displacement 
tensor, i.e. B l3k d k = e tJ . The last equation is convenient because it is identical 
to the tangent stiffness used in the local N-R for determining the slip, and 
the same code can be reused. The other equations can be further specified 
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dc 



1 

2 



§£% + 
df h i , 9fl 
as, b -r H 



dxi 



( 6 , 8 ) 



where C = <C» O'- 



Additionally, 



ch&\ 

~d> a 

<9<r 






da 



jk 



= n (X) n 



= (n 



iy 



( 6 . 9 ) 

( 6 . 10 ) 

( 6 . 11 ) 



The equations for the sliding-only case are as follows 



()c 






1 (df h , df 



d( s 



2 \ dx, 



In + 



dx 



-l, 



and 



— = (n ® l) s — sign (a) ji [n ® n) 



( 6 . 12 ) 



( 6 . 13 ) 



For the quantities corresponding to refer to Appendix A. 

The extra degrees of freedom may be condensed at the element level in the 
standard way, resulting in the final element stiffness matrix 
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6 . 1 . Enhancing Simulation Robustness 

In a comparative study of several different embedded discontinuity formula- 
tions, Jirasek [20] noted that even the most optimal enhanced finite element 
produces a global tangent stiffness matrix which is ill-conditioned and thus 
difficult to invert. Later, Jirasek and Belytschko in [21] compared the en- 
hanced finite element method with the well-known extended finite element 
method (XFEM). While both methods suffer from robustness issues arising 
from material softening, the XFEM has fewer drawbacks than the embedded 
discontinuity formulation. However, in the latter framework, because of the 
more straight-forward implementation and the fact that the added degrees of 
freedom within localized elements can be statically condensed at the element 
level, considerable effort has been made to develop methods to circumvent 
the aforementioned robustness issues. 

Although there are several remedies outlined in the literature, affording vary- 
ing levels of success, a particular method which seems to be especially promis- 
ing is an implicit-explicit integration technique, Impl-Ex, devised by Oliver 
and co-workers [28]. This method employs an implicit internal variable cal- 
culation at the end of the time step and on the subsequent time step, during 
the global N-R iteration, these values are used in a semi-implicit calculation 
of the stresses. The benefit of this method is that for the semi-implicit calcu- 
lation a positive definite stiffness matrix results which affords the simulation 
significant gains with regard to robustness. The conditioning of the stiffness 
matrix has also been shown to improve with the method [28]. 



Impl-Ex Integration Method 



For a given localized element, the internal variables such as plastic slip magni- 
tude are calculated implicitly at the end of the time step, after the convergent 
solution of the global displacements is obtained. Typically, on the subsequent 
step during the global N-R iteration, the stress in a localized element is cal- 
culated using, again, a fully implicit calculation. In contrast, for the Impl-Ex 
method, the stress is calculated semi-implicitly, which means that the pre- 
viously obtained implicit values from the prior time step are used as either 
semi-implicit or explicit values. The general form of the semi-implicit stress, 
which is recognizable as the standard predictor-corrector form used in stress 
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return plasticity algorithms, is given in Oliver by 



^n+l 



“I ~ C ■ As n ^_i AA n _|_iC • 



dg (<t w+ i) 
d&n+l 



(6.15) 



Using an explicit approximation of the plastic multiplier term we have 



A n +i — A n H — + AA n (6.16) 

L\I n 

where 

AA n = \ n — A„_i- (6-17) 

A n is the value of the previous time step’s plastic multiplier; it is calculated 
implicitly after convergence of the global N-R. Similarly, A n _i is the implicit 
value from the nth minus one time step. Thus, in the current time step that 
AA n is being used, it is an explicit value. 

Taking the derivative of the semi-implicit stress with respect to the current 
strain, we obtain the so-called effective algorithmic operator 



y'-feff 
°n+ 1 



da- 

ds 



n + 1 



■n+1 



= (I + AX n+1 c e : A 



n+l 



-1 



: c 



(6.18) 



A n+ \ 



d 2 g (<r n +i) 
dd n+1 (g) & n+1 



(6.19) 



(Refer to Appendix B for a derivation of the effective algorithmic operator.) 
In applying the Impl-Ex method to our particular formulation, the semi- 
implicit stress is calculated accordingly: 



<t „ +1 = <r„ + c" : Ae™' - c“ : (aC „ +1 ® V/'*) 



( 6 . 20 ) 



D. A. Weed et al./A Novel Post-localization Model 



25 



where 



Cn+i 






n + 



A*n±l ^ 

A t n v 



( 6 . 21 ) 



n+ 1 



For the effective operator calculation, we need the explicit form of the plastic 
potential term. If we define the direction of the slip on the deformation band 
as C 5 we have for the potential 



dg(& n + 1 ) 

<9dr n+1 



Cn+1 ® V/ n+1 



( 6 . 22 ) 



The explicit form is not necessary however, since £ n+1 <E> V/^ +1 is calculated 
during the previous time step. In this case, we freeze not only the magnitude 
but the direction of £ n+1 . I 11 the current time step, then, it is explicitly 
known and thus, A n+ \ = 0. Therefore, the jump in the displacement field is 
determined from the previous time step. Hence, the effective operator reduces 
to 



Cf +1 = I \ c e = c e (6.23) 

Hence, while there is a slight loss of accuracy in the added assumption, a 
simpler formulation is recovered that is linear in the case of linear elasticity. 
This fact tends to improve convergence rates in examples. On a conclud- 
ing note, while the Impl-Ex method provides tremendous gains in terms of 
numerical robustness and tractability, especially for problems with complex 
geometries, the drawback of this method is that accuracy and unconditional 
stability are sacrificed due to the semi-implicit calculation on the global step. 
Fortunately though, the remedy lies in simply choosing a sufficiently small 
time step such that the error is kept within acceptable bounds. In the numer- 
ical results section we show that with a small enough step size, the Impl-Ex 
produces results comparable to that of the fully implicit formulation. 
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7. Numerical Examples 

Four numerical examples are presented here as a comparison of the sliding 
model versus the combined sliding-opening model. 



7. 1 . Plane strain compression and tension of a column using a uniform 
load 

7.1.1. Plane strain compression test 

The first example is a pressure-confined column subjected to a uniform dis- 
placement (Figure 6). The loading is quasi-static; hence, inertial effects are 
neglected. The bottom of the specimen is fixed so that it cannot translate 
vertically. In addition, one node is fixed horizontally so that the specimen 
is allowed to expand laterally. A confining pressure of 25 MPa is applied on 
both sides of the specimen. The pressure is applied at the onset of the sim- 
ulation resulting in a uniform stress state throughout the mesh. This stress 
state causes all of the elements to localize in the same time step. This being 
the case, a seed element has to be chosen so that one dominant band will 
propagate throughout the mesh. When this initial element localizes, the band 
tracking algorithm will generate a single band from the seed element at the 
critical orientation and terminate at the opposite end of the mesh. 



Table 1 

Material parameters for plain strain compression sample 



Parameter 


Symbol 


Value 


Young’s Modulus 


E 


5500 MPa 


Poisson’s Ratio 


V 


0.25 


Cohesive Strength Parameter 


a 


8.034 MPa 


Friction Parameter 


p 


0.633 


Dilation Parameter 


b 


0.633 


Hardening Modulus 


H 


-10 MPa 


Localized Friction Coefficient 


g 


0.72 


Characteristic Slip Distance 


C 


0.5 mm 



The reaction force graph shows a characteristic elastic response as the mate- 
rial resists deformation, followed by a brief plastic response and then failure. 
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Fig 6: Geometry and loading for plane strain compression model. Simple 
constraints on the bottom allowing lateral expansion. 




(a) 




Fig 7: Localization lines and deformed shape of the plane. 






D. A. Weed et al./A Novel Post-localization Model 



28 




Fig 8: The graph shows the characteristic elastic and plastic responses as well 
as a stress drop and a constant frictional response. 



In quasi-brittle and rock-like materials, at low confining stresses, a narrow 
plastic region is expected due to the low ductility. After plasticity, failure 
occurs in which deformation in the localized elements is concentrated solely 
on the slip surface while the bulk material of the element unloads elastically. 
After the cohesion is fully degraded a constant friction force response is ob- 
served. Both the combined and sliding-only formulations produce identical 
results. This is to be expected since in a purely compressive state, both for- 
mulations simplify to the sliding Mohr-Coulomb law with friction. 



7.1.2. Tension test for plain- strain sample 

In this example, we reverse the direction of the displacement loading from 
the previous example and retain the confining forces on the specimen, thus 
creating a uniform stress in the vertical direction. As is expected, the re- 
versal of the displacement produces a slip line that is oriented 30°from the 
horizontal plane. From the reaction force profiles, the sliding-only formu- 
lation softens more quickly than the combined opening-sliding formulation 
due to the restricted kinematics of the system. A larger slip magnitude is 
required to obtain the same vertical displacement compared to the combined 
case, which opens in a more vertical pattern. When the cohesion degrades to 
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zero, both simulations fail to converge because the system becomes under- 
constrained (i.e. , the top half of the sample has no displacement constraints 
in the x-direction) . 




Fig 9: For the column in tension, the localization line is now oriented 30°from 
the horizontal plane. 



7.2. Sample with hole subjected to shearing load 

The next example, which features a slightly more complex geometrical figure, 
showcases the versatility that is gained from the additional opening degree of 
freedom. As mentioned previously, Foster et al [13] showed that the formula- 
tion involving the single degree of freedom led to spurious, mesh-dependent 
results. Namely, because the elements with different sliding directions are 
impeded by each other, stress build-up occurs around the localization zone. 

For the shearing example, in the case of the sliding-only formulation, multiple 
bands form (13b), which produce significant locking and subsequent harden- 
ing. Notice that due to this stress build-up, neither band is able to fully 
propagate. Figure 14 shows a distinctive kink in the reaction force which 
arises due to the hardening around the localization zone. In contrast, for 
the combined formulation one fully formed band (13a), on each side, occurs 
which is concomitant with a significant stress drop in the material. For higher 
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(c) opening-sliding 




dx (m) 




-2.0e-04 



-4.0e-04 



-5.10-04 



(b) sliding-only 




(d) sliding-only 



Fig 10: Comparison between the sliding-only and combined opening-sliding 
formulations for the pressure-confined column in tension. As expected, the 
sliding-only formulation shows more deformation in the horizontal direction 
whereas the opening-sliding model has greater deformation in the vertical 
direction. 
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Fig 11: Both simulations fail to converge once the cohesion reaches zero. The 
sliding-only case softens quicker because the cohesion only acts parallel to 
the deformation band whereas opening-sliding formulation presents cohesion 
both normal and parallel to the band. 



mesh refinements, after the cohesion completely degrades, there is a notice- 
able increase in the reaction force (16) which is due, not to geometric locking, 
but rather to more complex deformation patterns induced by the addition 
of more sampling points. This, in turn, produces an impediment to slippage, 
particularly at the corners where high compression is observed, resulting in 
the increase of shear resistance seen at the tail end of the softening curve. 
Nevertheless, a clear, single deformation band is observed in the deformation 
plot 15. Moreover, in the sliding-only formulation, the simulation is not able 
to complete the full loading schedule due to the interaction of multiple bands 
caused by the geometric locking impediment. The final case (13c) features el- 
ements which are embedded with a predefined strong discontinuity. This case 
was added to show that completely horizontal bands will produce a response 
in which the stress in the sample is completely released by the proliferation 
of fracture; in the absence of band curvature, the sample is allowed to freely 
translate in the horizontal direction. 

Figure 17 shows a comparison between the Impl-Ex and fully implicit inte- 
gration schemes for the shearing example with a hole, which employs only 
a sliding degree of freedom, comprised of 288 elements. For At = 1.0e-2 (a) 
the Impl-Ex shows noticeable oscillation, particularly at the peak and when 
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Table 2 

Material parameters for sample with hole 



Parameter 


Symbol 


Value 


Young’s Modulus 


E 


9000 MPa 


Poisson’s Ratio 


V 


0.15 


Cohesive Strength Parameter 


a 


8.034 MPa 


Friction Parameter 


p 


0.633 


Dilation Parameter 


b 


0.3165 


Localized Friction Coefficient 


l 1 


0.60 


Localized Hardening Modulus 


H 


0.0 MPa 


Characteristic Slip Distance 


c* 


0.5 mm 



d 




Fig 12: Boundary conditions for sample with a rectangular hole. There is a 
dominant shear displacement being applied to the top and a small vertical 
displacement which keeps the specimen from undergoing significant rotation 
after localization. 
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(a) sliding-opening 



(b) sliding-only 




(c) straight interface 

Fig 13: For the opening-sliding model (a) one dominant band is able to 
propagate whereas the sliding-only model (b) produces two bands both which 
experience locking. The bottom figure (c) displays the straight pre-embedded 
(interface) bands. 
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Fig 14: The sliding-only formulation shows locking which results in multiple 
bands, whereas the combined sliding-opening formulation shows a significant 
stress drop and a subsequent slight increase in stress due to some rotation of 
the sample at the end of the band Results are from a mesh comprised of 576 
elements. 




Fig 15: For the sliding-only formulation (a), two competing bands form and 
geometric locking ensues which produces a significant build up of bulk plas- 
ticity around the fracture zone. In the opening-sliding formulation (b), one 
dominant area of localization on each side of the specimen occurs, along 
which stress is released. 
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Fig 16: Convergence test for the plate with the hole run with the combined 
open-sliding formulation. The finest mesh shows some rise in the reaction 
force. This is not due to geometric locking but rather the compressive stresses 
at the corners of the hole. 
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the cohesion is significantly degraded. As the time step is refined, e.g. in the 
cases of At = 2.0e-3 (b) and At = 1.0e-3 (c) the Impl-Ex method is able to 
navigate these areas more adeptly. And in the final case of At = 2.0e-4 (d) 
the oscillation is mitigated to the point that it essentially becomes negligible. 
Though the time step needed for the Impl-Ex is significantly smaller than 
that of the fully implicit technique, it is far more robust. The fully implicit 
technique requires the use of elaborate time step-cutting schemes in order 
to reach convergence for more dense meshes, whereas the Impl-Ex does not. 
What is also notable is that the fully implicit is quite sensitive to the value of 
the critical slip parameter; if the material softens too quickly, this scheme has 
difficulty in capturing this quick descent and thus fails to converge at times. 
In contrast, because the Impl-Ex produces a positive-definite local stiffness 
matrix, it is not sensitive to how quickly the material softens. 



7.3. Plane strain slope stability 

This particular example is motivated by and thus similar to the one detailed 
in [32], This example was chosen because it presents a very complex mix of 
interactions, such as high compression, tension, friction and rotation. At the 
top of the slope is a rigid embankment which is comprised of elements that 
are only allowed to undergo elastic deformation; plasticity and localized be- 
havior are excluded. This is to prevent these elements from localizing so that 
softening behavior is solely relegated to the underlying material, which is the 
main area of interest. Figure 19 shows a comparison of the localization lines 
which illustrates the differing angle of descent between the two formulations. 
The combined formulation produces a more shallow angle. This can be under- 
stood in that the opening degree of freedom allows the material to undergo 
a significant amount of softening, over that of the single degree formulation, 
and hence less stress build-up occurs. In contrast, because geometric locking 
occurs in the sliding-only case, the resulting stress increase causes a steeper 
localization angle, more friction resistance and ultimately locking. 



In addition, a significant amount of rotation is seen in the combined case 
because the opening degree of freedom allows the material to move much 
more freely than in the sliding-only case. 



shear force, MN shear force, MN 
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Fig 17: Comparison between the Impl-Ex and the fnlly implicit formulations. 
Sliding-only formulation with 288 elements. 



Table 3 



Material parameters for slope 


stability 


example 


Parameter 


Symbol 


Value 


Young’s Modulus 


E 


10 MPa 


Poisson’s Ratio 


V 


0.4 


Cohesive Strength Parameter 


a 


40 KPa 


Friction Parameter 


p 


0.3 


Dilation Parameter 


b 


0.06 


Localized Friction Angle 


& 


10° 


Localized Hardening Modulus 


H 


0.0 MPa 


Characteristic Slip Distance 


c* 


0.4 m 
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20m 



Fig 18: Slope stability problem with point load applied on embankment. 




Fig 19: comparison of slip lines; the dotted black line is the opening-sliding 
case whereas the red line corresponds to the sliding-only case 



1000 



Sliding-only 




Fig 20: In the sliding-only case, there is significant friction resistance as well 
as locking and subsequent stress build-up. The opening-sliding formulation 
shows a large stress drop and resistance due to the slope rotation after com- 
plete cohesion degradation. 



8. Conclusions 

We have developed and implemented a constitutive model for localized soft- 
ening geomaterials. The model accounts for differing fracture energies in ten- 
sion and (mode II) shear, as well as friction along the interfaces. This model 
has been embedded in an enhanced strain finite element framework for mod- 
eling the propagation of discontinuities. 

In simulating the initiation and propagation of localized deformation, the 
combined opening-sliding constitutive model alleviates the spurious geomet- 
ric locking effect seen in the prior model which featured only a single slid- 
ing degree. Through several benchmark examples, the combined formulation 
displays expected softening as well as offers more accurate solutions, particu- 
larly in cases where tension is seen, even locally. Also, the formulation shows 
its versatility in large-scale problems as seen in the slope stability example. 
Lastly, implementation of the Impl-Ex technique affords the simulation code 
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d_y 

mm 9.04e-03 



-4.00e-01 



-5.67e-01 



(a) sliding-only 



d_y 

n 1.83e-01 
0.00e+00 



-4.00e-01 



-7.00e-01 




(b) opening-sliding 



Fig 21 
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(a) sliding-only 



d_x 

_ 3.406-02 

P 0.00e+00 



-4.00e-0l 



-5.20e-01 




(b) opening-sliding 



d_x 



■ 



4.50e-03 



-1.02e+00 



Fig 22 
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(b) opening-sliding 



Fig 23: Shear strain for a) sliding-only and b) opening-sliding cases. 
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significant improvements in convergence and robustness. With this in hand, 
the next logical step is to expand the simulations into the three-dimensional 
realm in order to simulate more realistic situations within a geomechanics 
context. 

Future studies will also include the combination of the this model with a 
plasticity more sophisticated model suited for geomaterials, which is detailed 
in [14, 24], The two models will be combined using the formulations outlined 
in [30]. In addition to this, the variable friction model in [13] can be incorpo- 
rated so as to more accurately capture the friction behavior of geomaterials. 
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Appendix A: Derivatives of the Yield Functions 

The gradients of the yield functions are as follows: 
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where 
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(A.7) 
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To determine the derivatives of the bulk stresses with respect to the jumps, 
recall 



o-n+i = o-n + c e : Ae conf - c e : (AC ® Wf h ) S 

Hence, 
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Appendix B: Derivation of the Effective Algorithmic Operator 



The effective algorithmic operator is defined as 

d&n+i 



/^ie ff 

L 'n+ 1 



de 



n + 1 



(B.l) 



by taking the derivative of the semi- implicit stress in eq. (6.15) with respect 
to the implicit strain we have the following: 



d&n+i 
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dcr n ^ _ d (e n+ i - e n ) 
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V da n+1 ) 



(B.2) 



where 



e de n _|_ 1 

C : h 

de n+ i 



AA n+ iC . A n _|_i . 



^ n +i 

de n+i 



A n+ i 
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Noting that c e : + 11+1 = c e : I = c e , we group terms and factor: 
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